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ABSTRACT 

We present a hydro dynamical simulation of the turbulent, magnetized, supernova (SN)-driven in¬ 
terstellar medium (ISM) in a stratified box that dynamically couples the injection and evolution of 
cosmic rays (CRs) and a self-consistent evolution of the chemical composition. CRs are treated as a 
relativistic fluid in the advection-diffusion approximation. The thermodynamic evolution of the gas 
is computed using a chemical network that follows the abundances of H+, H, H 2 , CO, C+, and free 
electrons and includes (self-)shielding of the gas and dust. We find that CRs perceptibly thicken the 
disk with the heights of 90% (70%) enclosed mass reaching > 1.5kpc {'> 0.2kpc). The simulations 
indicate that CRs alone can launch and sustain strong outflows of atomic and ionized gas with mass 
loading factors of order unity, even in solar neighborhood conditions and with a CR energy injection 
per SN of 10^^ erg, 10% of the fiducial thermal energy of an SN. The CR-driven outflows have moderate 
launching velocities close to the midplane (< lOOkms”^) and are denser (p ^ 10“^^ — 10“^^gcm“^), 
smoother, and colder than the (thermal) SN-driven winds. The simulations support the importance 
of CRs for setting the vertical structure of the disk as well as the driving of winds. 

Subject headings: cosmic rays — ISM: structure — ISM: jets and outflows — magnetohydrodynamics 
(MHD) — diffusion 


1. INTRODUGTION 


Galactic outflows are important for the dynamical 


and chemical evolution of galaxies (e.g. Veilleux, Gecil 


fc Bland-Hawthorn| 2005 Oppenheimer ST Dave 2006). 


I'hey can enrich galactic halos with metals and regulate 
the baryonic angular momentu m distribution in forming 
galaxies (e.g. Ubler et al. 2014|). Observations reveal that 


many actively star-forming galaxies i n the early Universe 


(z > 0.5) show wind signatures (e.g. 

Steidel et al. 12010; 

jNewman et al.|2012| Rubin et al. 12014 

[|. In the local Uni- 


(e.g. M82 or NGG 253). Milky Way (MW)-type galax¬ 
ies show galactic fountain features. Despite their impor¬ 
tance, the driving mechanisms of galactic winds are not 
precisely known. Galaxy-scale sub-grid modes including 
supernova (SN) explosions, radiation pressure, and stel¬ 
lar winds only partially succeed in explaining observed 
galactic wind properties. 

Recent numeri cal simulations of the galactic interstel- 
lar medium (ISM|Walch et al. 2015 Girichidis et al. 2015) 


have shown that the properties of SJN-driven outflows are 
tightly connected to details of the ISM structure and the 
environmental densities of SNe as one of the main driving 
mechanisms. Most previous calculations of the SN-driven 
ISM recover the net global quantities like mass fractions 
and volume filling factors of the differ ent thermal states 


of the gas close to the disk midplane (de Avillez & Bre 


itschwerdt||20051 IJoung & Mac Low||2006| Kim, Ostriker 

& Kim 2013t |Li et al.|2015p. However, the vertical dis- 

tribution of i 
perature stru 
numerical me 
Gosmic ray 

]he gas above the disk as well as the tem- 
.cture are not reproduced well by previous 

)dels (Henley et al.|2010 Wood et al.|2010). 

s (CRs) as a non-thermal component (see, 

e.g. Zweibel ! 

2013 Grenier, Black & Strong||2015P might 

change the structures in the ISIVI even more. Their en¬ 
ergy density is similar to the magnetic and thermal one 

(e.g., Gox||2005 Draine||2011|and references therein). As 

CRs do not ( 
provide a lor 
through the r 
cess relative 
large regions 
dient. Both 
dynamical dr 
From theor 
cut alone can 

200 I as etticiently as the thermal gas, they 
ig-lived energy reservoir. Their transport 
nedium can be described by a diffusion pro¬ 
to the gas, which allows CRs to populate 
of the disk creating a steady pressure gra- 
facts make CRs a fundamentally different 
iving mechanism in the ISM. 

‘etical considerations, a CR pressure gradi- 

generate outflows dBreitschwerdt, McKen- 

zie & Voelk|1993||Everett et al.|2008| Dorh & Breitschw- 


scale hydrodynamical simulations. Hanasz et al. (2013) 
neglect the thermal impact of SNe and show that CRs 
alone can launch fast winds from the magnetized ISM in 
high surface density disks (Egas = IOOMqPc”^). 


Booth 


et al. (|2Q13) assume isotropic diffusion without the mclu- 


sion of magnetic fields and compare MW and Small Mag- 
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ellanic Cloud (SMC)-type galaxies, finding that CRs no¬ 
ticeably increase the mass loading factor {r] = M/SFR) 
in SMC-like dwarf galaxies, but have a weaker impact in 
MW-type environmen ts. With a similar implementation, 
Salem & Bryan (2014) found stable outflows with 77 = 0.3 
but only for a high traction of CR energy injection per 
SN. The main conclusion from these studies is that CRs 
can diffuse out of the dense regions of star formation and 
generate stable vertical pressure gradients. This helps to 
lift even dense and cold gas above the disk midplane. 
However, these models do not resolve the ISM structures 
and phases. 

In this study, we present the first magnetohydro dy¬ 
namical numerical simulations of stratified boxes that 
dynamically include CRs and follow the chemical evolu¬ 
tion in the ISM. We aim for a more accurate description 
of the thermal and dynamical properties of a galactic disk 
with a focus on the vertical structure of the disk and the 
onset of outflows. 

2. NUMERICAL METHOD AND SIMULATION SETUP 

The simulations are performed with a modified version 
of the adaptive-mesh refinement code FLAS H 


m ver¬ 


sion 4 ( Fryxell et al.||2QQ'Q Dubey et al.||2QQ8 ). We solve 


HLL3R method (|Waagan|2009l iBouchut, Klingenberg & 

Waagan 

2010 Waagan, Lederrath Klingenberg||20Il|), 

extendec 

to a separate monoenergetic CR fluid similar 

to Yang 
presente* 

et al. (12012 

). Tests of this implementation are 

d inlCirichic 

Lis et al. (2014). 


merically is 


dpv 


| + v.M = o 


—+V^|pvv 


An 


de 


(e+Ptot)v- 


+ Vptot = pg 


B(B • v) 

Att 


pv • g + V • KVe^R + Mchem + ^inj 

— - V X (v X B) = 0 


de^ 


+ V • (CcrV) = 


-PcrV • V - 


dt 

V • (KVCc^r) + Qc 


( 1 ) 

( 2 ) 

(3) 

(4) 

(5) 


Here, p is the mass density, v is the velocity, and B is 
the magnetic field. The total energy density. 


e — pv^ ! 2 -h eth + 


(6) 


includes kinetic, thermal, CR, and magnetic contribu¬ 
tions. We evolve the CR energy density, separately. 
The total pressure is 

Ptot=Pth +PcR +Pmag (7) 

= (7 - l)eth + (7cr - l)ecR + (8) 

The closure relation for the system, the equation of state, 
combines the different contributions from CR and ther¬ 
mal pressure in an effective adiabatic index, 7 eff, 

T-Pth T 'Tcr-Pcr 


7eff = 


Pth +PcR 


(9) 


where we set 7 = 5/3 and 7 cr = 4/3 for gas and CRs, 
respectively. For the CR diffusion tensor, K, we as¬ 
sume a value of 10 ^^ cm^ s“^ parallel and 10 ^^ cm^ s~^ 
perpendicular to the magnetic field lines (e.g., Strong, 


Moskalenko & Ptuskin||2d67| |Nava & Gabici||20I^ 

An exact treatment of CR propagation and the cou¬ 
pling to the gas would include the generation of Alfven 
waves due to CR streaming and wave damping due to 
ion-neutral collisions and nonlinear Landau damping in- 
cluding anisotropic particle distribution functions (e.g. 
Kulsrud fc Pearc^ I969| [Ptuskin et ar]|I997[ [Breitschw- 


erdt, Dogiel fc V5ik||2QQ2| |Dorh fc Breitschwerdt||2QI2[ 

Uhlig et al.|2QI2 Wiener, Zweibel fc Oh|2Q13 ). However 


the detailed coupling of CRs to the gas - in particular, 
the weakly ionized gas - is an open question and beyond 
the scope of this Letter. For our model we assume that 
CRs are scattered sufficiently by small-scale background 
turbulence, which we do not resolve. We t herefore adopt 


Lerche 

(1985), see also 

Yang et al. ( 

2012); Hanasz et al. 


bution function. 

Radiative and chemical heating and cooling, 'Uchem, 
are computed using a chemical network that follows the 
abu ndances of H+, H, H 2 , C+, CO, and free electrons 
(see Glover & Clark 2QI2[ and references therein). We 
assume a constant CR ionization rate (3 x 10 “^^ s). We 
fol low (self-)shielding of the molec ular gas as described 
in Clark, Glover & Klessen| ( 2012 ). Photoionization is 
neglected. T'he thermal and CR energy injected by the 
SNe are included in Rinj and Qcr- T he external gravita- 
tional acceleration, g, is taken from |Kuijken fc Gilmore 


(1989). 

T'he stratified disk is initially atomic, at rest, and in 
pressure equilibrium with a scale hei ght of 60 pc and a 
gas surface density of I0 MqPc“^ (seejWalch et al.||20I5 
for further details) in a box of 2 x 2 x zli20kpc. Lor |z| < 
2.5 kpc, we adopt a cell size of 15.6 pc, which corresponds 
to an effective numerical resolution of 128 x 128 x 2560 
cells. For jzj > 2.5kpc, the adaptive-mesh refinement 
can coarsen the resolut ion. Motivated b y the Kennicutt- 
Schmidt (KS) relation ( Kennicutt|I998 ), we explode SNe 
at a fixed rate of 60 Myr~^k pc~^ assuming one massive 
star p er IOOM^t^. Similar to|de Avillez fc Breitschwerdt 


(2004) or Joung & Mac Low (12006), we assign 20% of 

J L „ _ 1„ I' UM -iT U, J Cno/ OAT X_TT 


the explosions to SIN type la and 80% to SN t ype H 
with scale heights of 325 pc and 50 pc, respectively (Tam- 
mann, Loefiler & Schroeder 1994). However, estimates 
of the displacement of runaway OB stars suggests that 
the effective scale height of type H SNe might be larger 
dLi et al.||20I5t. 

We separately vary the thermal energy injection per 
SN, Eth = /thF^SN, and the CR energy injection, Ecr = 
/cR^SN, with Esn = 10 ^^ erg. We discuss three differ¬ 
ent simulations: purely thermal SNe (/th = I.0,/cr = 
0.0, THl-CRO), purely CR SNe (/th = 0.0, /cr = O.I, 
THO-CRl), and the physically most realistic setup with 
both energy injections (/th = I.0,/cr = O.I, THl-CRl). 

3. MORPHOLOGICAL EVOLUTION 

In all simulations, the SNe induce the formation of fil¬ 
aments, clumps, and voids. The presence of CRs leads 
to a smoother morphological appearance. An overview of 
the disk structure after 250 Myr is depicted in Fig.[^with 
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Fig. 1.— Structure of the disk after 250 Myr of evolution. The left panel shows the run without CRs, the middle panel the run with 
purely CR feedback, and the right panel the simulation with both feedback energies. In the top and middle row, we plot the density edge-on 
and face-on in cuts through the center of the box. The bottom row depicts the CR energy density in the midplane. The disk in the run 
including CRs is significantly thicker compared to the purely thermal run. The thermal energy input of the SNe is primarily important in 
the denser part of the disk. CRs show their main impact in generating an extended atmosphere. The CR energy density is smooth in the 
midplane. Local variations are due to recent SNe and the resulting CR injection. Those variations disappear after less than a Myr due to 
diffusion. 
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Fig. 2. — Vertical profiles of the total gas density for all simula¬ 
tions. The arrows indicate the height of 90% enclosed mas s. A fit 
to the observed d ensity profile of the solar neighborhood (|Dickey| 
I&: Lockman||1990|) are shown in yellow. Thermal energy injection 
alone leads to a compact atomic gas distribution. Including CR 
feedback results in very extended distributions, which are much 
closer to the observed extent of the gas. The profiles indicate that 
CRs have their main impact at larger altitudes. 


cuts through the center of the box for the density edge- 
on (top) and face-on (middle) and the CR energy den¬ 
sity (bottom). From left to right, we present THl-CRO, 
THO-CRl, and THl-CRl. Both runs with CRs quickly push 
gas away from the midplane and form an extended atmo¬ 
sphere, whereas in the simulation with only thermal SN 
feedback, the disk remains compact without a noticeable 
extended layer of gas above the disk. The CR energy 
density shows relatively variations of about one order of 
magnitude. High CR energy peaks correspond to recent 
SNe and disappear quickly due to diffusion. 

In Fig. we show the vertical profile of the density 
for all simulations at t = 250 Myr. The arrows indicate 
the heights of 90% enclosed mass. The inclusion of CRs 
increases i^ 9 o% from ^ 30 pc to ^ 1500 pc with similar 
values for both runs including CRs. The additional ther¬ 
mal energy injection from SNe in run THl-CRl changes 
the profile up to a height of ^ 300 pc. At larger alti¬ 
tudes, the profiles for THO-CRl and THl-CRl are basically 
indistinguishable. For compar ison, we include the verti- 
cal density profile of the MW (Dickey & Lockman||1990) 
showing reasonable agreement with both CR runs. 

The time evolution of the disk thickness is shown in 
Fig. [^including 70% (solid line) and 90% (dashed line) of 
enclosed mass. Run THl-CRO shows little variation of the 
disk over the entire simulation time. For the other runs, 
the CR energy density can push a significant fraction of 
the gas to heights of up to ^ 400 pc in the beginning of 
the simulation before the initially ordered magnetic field 
becomes tangled (Parker instability) and the CR energy 
quickly diffuses along the opened up field lines. Opening 
magnetic field lines enable the reduction of CR pressure 
via the escape of excessive CR gas initially trapped in the 
horizontal magnetic field. The bulk of the disk thickens 
slowly over time whereas the envelope expands 

quickly 


4. OUTFLOWS 

The properties of the outflows are plotted in Fig. 
with the total outflow rate (top), the ratio of atomic to 
ionized hydrogen (middle), and the temperature (bot¬ 
tom), averaged at Tlkpc. Included is all the gas with 


velocities pointing away from the midplane. All simu¬ 
lations drive outflows at a roughly constant rate after 
^ 50 Myr. Run THl-CRO reaches a mass loading factor, 
T] = M/SFR ^ 0.1 — 0.2 at Ikpc above the midplane. 
Including CRs increases r] to order unity. This indicates 
that in our setups the main driver for outflows is CRs 
rather than the thermal contribution from the SNe. The 
chemical composition, however, is affected by the ther¬ 
mal contribution. The outflowing gas in the simulation 
with purely thermal SN feedback is very hot (10^ K; bot¬ 
tom panel of Fig.|^ and composed entirely of ionized hy¬ 
drogen H+ (not visible). In the other two cases with CRs, 
the outflowing gas is warm (^ 10^ K) and mainly com¬ 
posed of atomic hydrogen with Mh/^h+ ^5 — 10. The 
inclusion of photoionization feedback is likely to decrease 
this value, but will not significantly alter the tempera¬ 
ture of the outflowing gas or the very low volume-filling 
fractions of the hot gas. As soon as the CR pressure 
gradient has built up, significant fractions of neutral and 
warm gas can be accelerated away from the disk mid¬ 
plane (see, e.g. jSalem & Bryan||2014|). 

The connection between density and outflow velocity 
is depicted in Fig. For THl-CRO (top panel), the SNe 
explode partially in regions of very low density, can thus 
expand to large radii, and push gas out of the midplane 
at very high velocities of up to ^ 500kms“^ at very low 
densities below p ^ 10“^^gcm“^. If CRs are included, 
gas is slowly transported away from the dense regions 
in the midplane to form an extended atmosphere, whose 
flow of gas becomes perceptibly denser and overall slower 
compared to THl-CRO. At the low heights investigated 
here, the gas does not reach escape velocity. Again, this 
is because the CRs quickly diffuse throughout the volume 
and build up a stable, slowly varying energy (pressure) 
gradient along the 2 ;-direction that slowly but continu¬ 
ously accelerates the gas, in contrast to the short-lived 
locally acting thermal energy input. This underlines 
the fundamental difference of galactic winds driven by 
SNe and CRs. The comparison between THO-CRl (mid¬ 
dle) and the combined energy input (THl-CRl) shows 
only little impact on the bulk density of the outflow 
(p ^ 10“^^ — 10“^^gcm“^). The contribution of the 
thermal component in the SN driving increases the peak 
velocities slightly from ^ 50kms~^ to ^ lOOkms”^. 
However, the dense atmosphere around the midplane, 
in which the SNe explode, reduces their net impact on 
the velociti es compared to THl-CRO (see also |Girichidis 
et al.|[2M5 ). For observationa l signa tures in tne soft X- 


ray emission, see 


Peters et al. (2015) 


5. DISCUSSION 

We find that CRs are able to perceptibly increase 
the thickness of the disk and launch and sustain strong 
outflows with mass loading factors of order unity. For 
the simulations presented here, the CR-driven winds are 
slower (50 — lOOkrns”^) and denser than the (pure ther¬ 
mally) SN-driven winds. This is qualit a tively consistent 
with t he results by Booth et al. (2013|). Salem & Bryan 
(|2014|) report mass loading factors below unity tor their 
hducial model with a total CR fraction p er SN of 30% 
rather than our 10%. The simulations of lHanasz et al.l 
(2013) also show mass loading factors above unity but 
for massive galactic disks with high SN rates. We note. 
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Fig. 3. — Time evolution of the vertical gas distribution. Color coded is the total gas fraction. Overplotted are the heights of 90% 
dashed line) and 70% {Hyo%, solid line) enclosed total mass. In the purely thermal run (left) the disk does not show variations over time. 
Including CRs leads to a slow expansion of bulk of the disk and a fast expansion of the envelope 
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Fig. 4. — Outflow rate (top), the ratio of atomic to ionized hy¬ 
drogen (middle), and the temperature (bottom) at 2 ; = ±1 kpc. All 
simulations show outflows. The purely thermal SNe have outflow 
rates lower than the star formation rate. Including CRs, the mass 
loading factor increases to unity. The outflows driven by purely 
thermal SNe are hot (~ 10^ K) and consist of ionized hydrogen. 
The CR-driven counterparts are colder (10"^ K) and mainly com¬ 
posed of atomic hydrogen. 
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Fig. 5. — Histograms of the outflows {vzZ > 0) as a function of 
density and velocity. Purely thermal SN feedback (top) leads to 
fast low-density outflows. Including CRs (middle and bottom) the 
transported gas is two orders of magnitude denser and a factor of 
a few slower with broader distributions 
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however, that an important parameter of the model is the 
exact value of the assumed CR diffusion coefficient. Its 
value con trols both the terminal wind speed and mass- 
loss rate ( Dorfi fc Breitschwerdt|[2M2 ). For small values 
of K, mass-loss rates are higher at lower speeds, whereas 
large diffusion coefficients result in lower mass loading 
factors but higher velocities. We found a similar behav¬ 
ior in tests with varying diffusion tensor for the setup 
described here. Full details will be provided in a follow¬ 
up paper. 

Previous models of stratified disks have shown that 
magnetic fields he lp to increase the thickness of the disk 
(Hill et al. 2012), but that they alone are not strong 
enough to reach observed scale heights. The combined 
MHD-CR simulation presented here seems a promising 
solution to this problem. 

We also note that the impact of purely thermal SNe 
strongly depends on the positioning of the SNe and there¬ 
fore might be resolution dependent. In particular, the 
SNe that explode in density peaks have a weaker net im¬ 


pact (Gatto et al. 2015 Walch et al. 2015 Girichidis et al. 


2015), and thus the outhow rates in the purely thermal 


run, THl-CRO, might be lower limits. This weak SN effi¬ 
ciency might also be responsible for small fractions of hot 
gas in the CR runs. However, we find that the inclusion 
of CRs results in outflows that are compatible with the 
observations. This suggests that CRs will always play an 
important role in the matter cycle of galaxies. 


6. CONCLUSIONS 

We present the first ISM simulations that dynamically 
couple CRs with the chemodynamical evolution of the 
magnetized ISM. The CRs are included in an advection- 
diffusion approach with an anisotropic diffusion tensor. 
We follow the formation of density structures, the verti¬ 
cal distribution of gas in a galactic disk, and the launch¬ 
ing of outflows in three different runs varying the SN 
feedback between purely thermal, purely CR, and com¬ 
bined thermal and CR feedback. Our conclusions can be 
summarized as follows: 


1 . Including CRs thickens the galactic disk. The 
height of 90% enclosed total mass is found to be 
^ 1.5 kpc in the case of 10% CR energy injection 
per SN after 250 Myr and to increase continuously. 
Comparison with the vertical density distribution 
in the MW indicates good agreement. 


2 . We find that CRs quickly lead to the formation 


of a warm and neutral galactic atmosphere pro¬ 
viding a mass reservoir for galactic winds and out¬ 
flows. Whereas the thermal contribution of the SNe 
mainly shapes the disk close to the midplane, the 
additional CR energy shows the strongest impact 
above the disk and in the halo. 

3. All simulations drive gas out of the midplane 
with little variation over time. For purely ther¬ 
mal SN feedback, the outflows are hot and com¬ 
posed of mainly ionized hydrogen with rates be¬ 
low the star formation rate. They are fast (up 
to ^ a few lOOkms”^) with low densities (p < 
10 “^^gcm“^). CRs alone can drive outflows with 
mass loading factors of order unity, which are warm 
(10^ K) and mainly composed of atomic hydrogen. 
They are a factor of a few slower 10 — 50 km s~ ) 

and 1 — 2 orders of magnitude denser (p ^ 10“^^ — 
10-25 gcm“3) compare to their thermally driven 
counterparts. 
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